{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению. Сессия № 3\n",
    "\n",
    "### <center> Автор материала: Михаил Каменев "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Индивидуальный проект по анализу данных </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**План исследования**\n",
    " - Описание набора данных и признаков\n",
    " - Первичный анализ признаков\n",
    " - Первичный визуальный анализ признаков\n",
    " - Закономерности, \"инсайты\", особенности данных\n",
    " - Предобработка данных\n",
    " - Создание новых признаков и описание этого процесса\n",
    " - Кросс-валидация, подбор параметров\n",
    " - Построение кривых валидации и обучения \n",
    " - Прогноз для тестовой или отложенной выборки\n",
    " - Оценка модели с описанием выбранной метрики\n",
    " - Выводы\n",
    " \n",
    " Более детальное описание [тут](https://goo.gl/cJbw7V)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "import scipy\n",
    "from statsmodels.stats.weightstats import *\n",
    "\n",
    "from sklearn.linear_model import RidgeCV, Ridge, Lasso, LassoCV\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import train_test_split, KFold\n",
    "from sklearn.preprocessing import StandardScaler, LabelBinarizer, PolynomialFeatures\n",
    "from sklearn.metrics import mean_absolute_error, make_scorer\n",
    "from xgboost import XGBClassifier\n",
    "from hyperopt import fmin,tpe, hp, STATUS_OK, Trials\n",
    "import xgboost as xgb\n",
    "from sklearn.model_selection import learning_curve, validation_curve\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 1. Описание набора данных и признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Датасет содержит информацию о 53940 бриллиантах. По некоторым характеристикам (об этом позже) будем предсказывать стоимость. Данные можно скачать <a href='https://www.kaggle.com/shivam2503/diamonds/data'>здесь</a>.\n",
    "\n",
    "С точки зрения бизнеса, ценность задачи понятна - по характеристикам бриллианта предсказать, сколько долларов за него можно получить. От бизнеса я далёк, поэтому интерес чисто спортивный: разобраться, какие характеристики и как влияют на стоимость этих камешков =)\n",
    "\n",
    "<b>Признаки</b>\n",
    "- carat - вес бриллианта в каратах, вещественный\n",
    "- cut - качество огранки, категориальный. Принимает пять возможных значений: Fair, Good, Very Good, Premium, Ideal\n",
    "- color - \"цвет\" бриллианта. Категориальный признак, принимает значения J,I,H,G,F,E,D (от худшего (J) к лучшему (D))\n",
    "- clarity - чистота бриллианта. Категориальный признак, принимает значения I1 (худший), SI2, SI1, VS2, VS1, VVS2, VVS1, IF (лучший)\n",
    "- x,y,z - три признака, характеризущие размеры бриллианта, вещественные\n",
    "- depth - признак, который высчитывается на основе трех предыдущих по формуле 2 * z / (x + y), вещественный\n",
    "- table - отношение ширины верхней грани бриллианты к его максимальной ширине, в процентах\n",
    "\n",
    "\n",
    "<b>Целевой признак</b>: price - стоимость бриллианта в долларах\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 2. Первичный анализ признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#загружаем dataset\n",
    "diamonds_df = pd.read_csv('../../data/diamonds.csv')\n",
    "diamonds_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diamonds_df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что масштабы признаков отличаются. В дальнейшем нужно будет применить StandartScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diamonds_df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В данных отсутствуют пропуски. Итого, имеется 6 вещественных, 1 целочисленный (unnamed: 0 не считаем) и 3 категориальных признака."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Анализ целочисленных и вещественных признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "real_features = ['carat', 'depth', 'table', 'x', 'y', 'z','price']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Изучим корреляцию вещественных признаков и целевой переменной\n",
    "sns.heatmap(diamonds_df[real_features].corr(method='spearman'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Признаки carat, x,y,z имеют большую корреляцию, как между собой, так и с целевой переменной, что не удивительно. При этом, корреляция целевой переменной и признаков depth, table почти отсутствует"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Анализ категориальных признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat_features = ['cut','color','clarity']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 10))\n",
    "\n",
    "for idx, feature in enumerate(cat_features):\n",
    "    sns.countplot(diamonds_df[feature], ax=axes[idx % 3], label=feature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Реальные значения категориальных признаков не отличаются от тех, что заявлены в описании. Кроме того, видно, что уникальных значений не много, так что One Hot encoding должен отлично сработать. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Анализ целевого признака"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(diamonds_df['price'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Распределение имеет тяжелый правый хвост. Применим логарифмирование."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(diamonds_df['price'].map(np.log1p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Помогло это не сильно: получилось бимодальное распределение. Зато хвост исчез =) Для наглядности, построим QQ график"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.probplot(diamonds_df['price'], dist=\"norm\", plot=plt);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Выводы\n",
    "- Вещественные признаки (carat, depth, table, x, y, z) масштабируем\n",
    "- К категориальным признакам ('cut','color','clarity') применяем one hot encoding\n",
    "- Целевую переменную логарифмируем"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 3. Первичный визуальный анализ признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Анализ целочисленных и вещественных признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Начнем с построения гистограмм вещественных признаков\n",
    "fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))\n",
    "\n",
    "for idx, feature in enumerate(real_features[:-1]): #price рисовать не будем\n",
    "    sns.distplot(diamonds_df[feature], ax=axes[idx // 3, idx % 3], label=feature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Распределение признаков depth, table, y, z отдаленно, но напоминает колокол. У depth хвосты тяжеловаты для нормального распределения; carat и table скорее бимодальные. Кроме того, у них тяжелые правые хвосты, так что np.log1p не помешает. По графикам выше не видно выбросов. Проверим, что это действительно так, с помощью boxplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))\n",
    "\n",
    "for idx, feature in enumerate(real_features[:-1]): #price рисовать не будем\n",
    "    sns.boxplot(diamonds_df[feature], ax=axes[idx // 3, idx % 3], orient='v')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Каких-либо серьезных аномалий в рассматриваемых данных нет. На всякий случай посмотрим бриллиант с y=60, z = 32 и carat > 4. Если он стоит дорого, то за выброс его считать не будем."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diamonds_df[diamonds_df['y'] > 55].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diamonds_df[diamonds_df['z'] > 30].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diamonds_df[diamonds_df['carat'] > 4].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что это просто очень дорогие камни. Посмотрим, как рассматриваемые признаки взаимосвязаны с целевой перменной"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.pairplot(diamonds_df[real_features], diag_kind=\"kde\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- вес бриллианта показывает степенную зависимость от его размеров\n",
    "- depth и table почти никак не взаимосвязаны с остальными признаками, в том числе и целевым\n",
    "- x,y,z связаны между собой линейно\n",
    "- цена линейно зависит от размеров\n",
    "- зависимость между ценой и весом сложно назвать линейной, но монотонный тренд есть"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Анализ категориальных признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим, как целевая переменная зависит от категориальных признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# цвет бриллианта\n",
    "fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(16, 10))\n",
    "\n",
    "for idx, (color, sub_df) in  enumerate(pd.groupby(diamonds_df, 'color')): \n",
    "    ax = sns.distplot(sub_df['price'], kde=False,  ax=axes[idx // 3, idx % 3])\n",
    "    ax.set(xlabel=color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Распределения для всех значений цветов имеют тяжелый правый хвост и не сильно отличаются друг от друга."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# чистота бриллианта\n",
    "fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(16, 10))\n",
    "\n",
    "for idx, (clarity, sub_df) in  enumerate(pd.groupby(diamonds_df, 'clarity')): \n",
    "    ax = sns.distplot(sub_df['price'], kde=False,  ax=axes[idx // 3, idx % 3])\n",
    "    ax.set(xlabel=clarity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Хвосты у всех тяжелые, но у SI1,SI2 присутствуют дополнительные пики в районе 5000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# качество огранки\n",
    "fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))\n",
    "\n",
    "for idx, (cut, sub_df) in  enumerate(pd.groupby(diamonds_df, 'cut')): \n",
    "    ax = sns.distplot(sub_df['price'], kde=False,  ax=axes[idx // 3, idx % 3])\n",
    "    ax.set(xlabel=cut)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И снова пики в районе 5000 (у Good и Premium). А в целом графики похожи."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Нарисуем boxplot для каждого значения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(16, 10))\n",
    "\n",
    "# Отобразим строки в числа в порядке от худшего к лучшему. Так удобнее на графике смотреть\n",
    "df = diamonds_df.copy()\n",
    "df['color'] = df['color'].map({'J': 0, 'I': 1, 'H': 2, 'G': 3, 'F': 4, 'E': 5, 'D': 6})\n",
    "df['clarity'] = df['clarity'].map({'I1': 0, 'SI2': 1, 'SI1': 2, 'VS2': 3, 'VS1': 4, 'VVS2': 5, 'VVS1': 6, 'IF': 7 })\n",
    "df['cut'] = df['cut'].map({'Fair': 0, 'Good': 1, 'Very Good': 2, 'Premium': 3, 'Ideal': 4})\n",
    "\n",
    "for idx, feature in enumerate(cat_features):\n",
    "    sns.boxplot(x=feature, y='price',data=df,hue=feature,  ax=axes[idx])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Тут уже интереснее. Начнем с огранки. Видно, что медиана максимальна для Very Good и Premium. Для ideal медианное значение цены гораздо меньше. Аналогичные наблюдения можно сделать для цвета и чистоты. Возможно, бриллианты с наилучшими свойствами на очень большие, и, соответсвенно, их цена ниже. Проверим это."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(16, 10))\n",
    "\n",
    "for idx, feature in enumerate(cat_features):\n",
    "    sns.boxplot(x=feature, y='carat',data=df,hue=feature,  ax=axes[idx])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Действительно, медианное значение веса для бриллиантов с очень хорошими характеристиками меньше, чем для бриллиантов с плохими харакетристиками. Напоследок, посмотрим сколько бриллиантов с той или иной харакеристикой присутствует в данных."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 10))\n",
    "\n",
    "for idx, feature in enumerate(cat_features):\n",
    "    sns.countplot(df[feature], ax=axes[idx % 3], label=feature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что очень мало камней с плохой огранкой. Также мало камней с плохими цветовыми харакеристиками. Но и не очень много с идеальными. Распределение чистоты камня напоминает лапласовское распределение."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 4. Закономерности, \"инсайты\", особенности данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Основные выводы по предыдущим пунктам:\n",
    "- Ключевые признаки для прогнозирования: вес и размеры бриллианта (carat, x, y, z). По графикам видно, что есть монотонная зависимость этих признаков и цены. Что логично\n",
    "- Признаки depth и table почти не влияют на стоимость камня\n",
    "- Исключительно по категориальным признакам сложно что-либо сказать о целевой переменной. Однако видно, что чем лучше бриллиант с точки зрения этих признаков, тем больше вероятность того, что он будет не очень большого размера\n",
    "- Выбросы в данных отсутствуют\n",
    "- Так как у целевой переменной очень тяжелый правый хвост, в качестве метрики будем использовать среднюю абсолютную ошибку, а не квадратичную.\n",
    "- Видно, что зависимость от ключевых признаков близка к линейной. Поэтому в качестве бейзлайна будем использовать линейную регрессию.\n",
    "- Более того, признаков не так уж и много, поэтому будем рассматривать также случайный лес и градиентный бустинг (тут он должен затащить =)). А случайный лес инетересен исключительно для сравнения с бустингом"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 5. Предобработка данных "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Для начала, выделим выборку для тестирования\n",
    "X = diamonds_df.drop(['price'], axis=1).values[:,1:] # отсекаем индекс\n",
    "y = diamonds_df['price']\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.3,random_state=4444, shuffle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# признаки с индексами 1,2,3 категориальные. Применим к ним ohe\n",
    "label_bin = LabelBinarizer()\n",
    "X_train_cut_ohe = label_bin.fit_transform(X_train[:,1])\n",
    "X_test_cut_ohe = label_bin.transform(X_test[:,1])\n",
    "                            \n",
    "X_train_color_ohe = label_bin.fit_transform(X_train[:,2])\n",
    "X_test_color_ohe = label_bin.transform(X_test[:,2])\n",
    "                                                    \n",
    "X_train_clarity_ohe = label_bin.fit_transform(X_train[:,3])\n",
    "X_test_clarity_ohe = label_bin.transform(X_test[:,3])  \n",
    "\n",
    "# carat, x и целевую переменную логарифмируем\n",
    "log_vect = np.vectorize(np.log1p)\n",
    "X_train_сarat_log = log_vect(X_train[:,0]).reshape(-1,1)\n",
    "X_test_сarat_log  = log_vect(X_test[:,0]).reshape(-1,1)\n",
    "X_train_x_log = log_vect(X_train[:,6]).reshape(-1,1)\n",
    "X_test_x_log  = log_vect(X_test[:,6]).reshape(-1,1)\n",
    "y_train_log = log_vect(y_train)\n",
    "y_test_log = log_vect(y_test)\n",
    "\n",
    "# масштабириуем вещественные признаки\n",
    "scaler = StandardScaler()\n",
    "X_train_real = np.hstack((X_train_сarat_log, X_train_x_log, X_train[:,[7,8,4,5]]))\n",
    "X_test_real = np.hstack((X_test_сarat_log, X_test_x_log, X_test[:,[7,8,4,5]]))\n",
    "X_train_real_scaled = scaler.fit_transform(X_train_real)\n",
    "X_test_real_scaled = scaler.transform(X_test_real)\n",
    "\n",
    "# В качестве дополнительных признаков будем рассматривать полиномиальные признаки\n",
    "#Данные признаки должны улучшить качество линейной модели.\n",
    "\n",
    "X_train_additional = PolynomialFeatures().fit_transform(X_train_real)\n",
    "X_test_additional = PolynomialFeatures().fit_transform(X_test_real)\n",
    "X_train_additional_scaled = scaler.fit_transform(X_train_additional)\n",
    "X_test_additional_scaled = scaler.transform(X_test_additional)\n",
    "\n",
    "\n",
    "\n",
    "# Объединяем все преобразованные признаки\n",
    "X_train_transformed = np.hstack((X_train_real_scaled,X_train_cut_ohe, X_train_color_ohe, X_train_clarity_ohe))\n",
    "X_test_transformed = np.hstack((X_test_real_scaled,X_test_cut_ohe, X_test_color_ohe, X_test_clarity_ohe))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 6. Создание новых признаков и описание этого процесса"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Смотри предыдущий пункт"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 7. Кросс-валидация, подбор параметров"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Рассмотрим сначала линейную модель. Данные разделим на 5 фолдов. С помощью RidgeCV и LassoCV будем оптимизировать силу регуляризации."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# функция потерь для рассматриваемой задачи. Ошибку смотрим на исходных данных\n",
    "def mean_absolute_exp_error(model, X,y):\n",
    "    return -mean_absolute_error(np.expm1(model.predict(X)), np.expm1(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv = KFold(n_splits=5, shuffle=True, random_state=4444)\n",
    "alphas = np.logspace(-5,2,100)\n",
    "ridge_cv = RidgeCV(alphas=alphas, scoring=mean_absolute_exp_error, cv=cv)\n",
    "lasso_cv = LassoCV(alphas=alphas, cv=cv, random_state=4444)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ridge_cv.fit(X_train_transformed, y_train_log)\n",
    "lasso_cv.fit(X_train_transformed, y_train_log)\n",
    "print('Optimized alpha: Ridge = %f, Lasso = %f' % (ridge_cv.alpha_, lasso_cv.alpha_))\n",
    "score_ridge = mean_absolute_error(y_test, np.expm1(ridge_cv.predict(X_test_transformed)))\n",
    "score_lasso = mean_absolute_error(y_test, np.expm1(lasso_cv.predict(X_test_transformed)))\n",
    "print('Ridge regression score = %f' % score_ridge)\n",
    "print('Lasso regression score = %f' % score_lasso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Оба метода показали схожий результат. Что будет, если мы добавим новые признаки?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_transformed_add = np.hstack((X_train_transformed, X_train_additional_scaled))\n",
    "X_test_transformed_add = np.hstack((X_test_transformed, X_test_additional_scaled))\n",
    "ridge_cv.fit(X_train_transformed_add, y_train_log)\n",
    "lasso_cv.fit(X_train_transformed_add, y_train_log)\n",
    "print('Optimized alpha: Ridge = %f, Lasso = %f' % (ridge_cv.alpha_, lasso_cv.alpha_))\n",
    "score_ridge = mean_absolute_error(y_test, np.expm1(ridge_cv.predict(X_test_transformed_add)))\n",
    "score_lasso = mean_absolute_error(y_test, np.expm1(lasso_cv.predict(X_test_transformed_add)))\n",
    "print('Ridge regression score = %f' % score_ridge)\n",
    "print('Lasso regression score = %f' % score_lasso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ошибка значительно уменьшилась. Построим кривые валидации и обучения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# код из статьи на хабре\n",
    "\n",
    "\n",
    "def plot_with_err(x, data, **kwargs):\n",
    "    mu, std = data.mean(1), data.std(1)\n",
    "    lines = plt.plot(x, mu, '-', **kwargs)\n",
    "    plt.fill_between(x, mu - std, mu + std, edgecolor='none',\n",
    "    facecolor=lines[0].get_color(), alpha=0.2)\n",
    "\n",
    "model = Ridge(random_state=4444) \n",
    "alphas = np.logspace(1,2,10) + 10 # Если коэффициент регуляризации мал, то значения получаются заоблочными\n",
    "val_train, val_test = validation_curve(model, X_train_transformed_add, y_train_log,'alpha', alphas, cv=cv,scoring=mean_absolute_exp_error)\n",
    "plot_with_err(alphas, -val_train, label='training scores')\n",
    "plot_with_err(alphas, -val_test, label='validation scores')\n",
    "plt.xlabel(r'$\\alpha$'); plt.ylabel('MAE')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Судя по кривым валидации, модель недообучилась: ошибки лежат близко друг к другу. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# код из статьи на хабре\n",
    "    \n",
    "def plot_learning_curve(model, X,y):\n",
    "    train_sizes = np.linspace(0.05, 1, 20)\n",
    "      \n",
    "    N_train, val_train, val_test = learning_curve(model,X, y, train_sizes=train_sizes, cv=5,scoring=mean_absolute_exp_error, random_state=4444)\n",
    "    plot_with_err(N_train, -val_train, label='training scores')\n",
    "    plot_with_err(N_train, -val_test, label='validation scores')\n",
    "    plt.xlabel('Training Set Size'); plt.ylabel('MAE')\n",
    "    plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Ridge(alpha=52.140083,random_state=4444)\n",
    "plot_learning_curve(model, X_train_transformed_add, y_train_log)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Кривые лежат близко друг к другу почти с самого начала. Вывод: наблюдений у нас достаточно, нужно двигаться в сторону усложнения модели"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Случайный лес"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Случайный лес должен хорошо работать \"из коробки\". Поэтому будем оптимизировать только число деревьев.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "model = RandomForestRegressor(n_estimators=100, random_state=4444)\n",
    "n_estimators = [10,25,50,100,250,500,1000]\n",
    "val_train, val_test = validation_curve(model, X_train_transformed, y_train_log,'n_estimators', n_estimators, cv=cv,scoring=mean_absolute_exp_error)\n",
    "    \n",
    "plot_with_err(n_estimators, -val_train, label='training scores')\n",
    "plot_with_err(n_estimators, -val_test, label='validation scores')\n",
    "plt.xlabel('n_estimators'); plt.ylabel('MAE')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что начиная с 200 деревьев качество практически не изменяется. Поэтому в качестве еще одной модели будем рассматривать случайный лес именно с таким количеством деревьев."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forest_model = RandomForestRegressor(n_estimators=200, random_state=4444)\n",
    "forest_model.fit(X_train_transformed, y_train_log)\n",
    "forest_prediction = np.expm1(forest_model.predict(X_test_transformed))\n",
    "score = mean_absolute_error(y_test, forest_prediction)\n",
    "print('Random forest score: %f' % score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# посмотрим на важность признаков\n",
    "np.argsort(forest_model.feature_importances_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Первые четыре столбца обучающей выборки соответствуют признакам carat, x,y,z. Как и предполагалось в начале, 3 из 4 этих признаков имеют наибольшую важность для модели"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Построим, также, кривую обучения\n",
    "plot_learning_curve(model, X_train_transformed, y_train_log)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "График выходит на полку, так что больше данных нам не нужно"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### boosting. А что boosting?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_boosting, X_valid_boosting, y_train_boosting, y_valid_boosting = train_test_split(\n",
    "    X_train_transformed, y_train_log, test_size=0.3, random_state=4444)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def score(params):\n",
    "    from sklearn.metrics import log_loss\n",
    "    print(\"Training with params:\")\n",
    "    print(params)\n",
    "    params['max_depth'] = int(params['max_depth'])\n",
    "    dtrain = xgb.DMatrix(X_train_boosting, label=y_train_boosting)\n",
    "    dvalid = xgb.DMatrix(X_valid_boosting, label=y_valid_boosting)\n",
    "    model = xgb.train(params, dtrain, params['num_round'])\n",
    "    predictions = model.predict(dvalid).reshape((X_valid_boosting.shape[0], 1))\n",
    "    score = mean_absolute_error(np.expm1(y_valid_boosting), np.expm1(predictions))\n",
    "   # score = mean_absolute_error(y_valid_boosting, predictions)\n",
    "    print(\"\\tScore {0}\\n\\n\".format(score))\n",
    "    return {'loss': score, 'status': STATUS_OK}\n",
    "\n",
    "def optimize(trials):\n",
    "    space = {\n",
    "             'num_round': 200,\n",
    "             'learning_rate': hp.quniform('eta', 0.05, 0.5, 0.005),\n",
    "             'max_depth': hp.quniform('max_depth', 3, 14, 1),\n",
    "             'min_child_weight': hp.quniform('min_child_weight', 1, 10, 1),\n",
    "             'subsample': hp.quniform('subsample', 0.5, 1, 0.05),\n",
    "             'gamma': hp.quniform('gamma', 0.5, 1, 0.01),\n",
    "             'colsample_bytree': hp.quniform('colsample_bytree', 0.4, 1, 0.05),\n",
    "             'eval_metric': 'mae',\n",
    "             'objective': 'reg:linear',\n",
    "             'nthread' : 4,\n",
    "             'silent' : 1,\n",
    "             'seed': 4444\n",
    "             }\n",
    "    \n",
    "    best = fmin(score, space, algo=tpe.suggest,trials=trials, max_evals=100)\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Оптимизация параметров\n",
    "trials = Trials()\n",
    "best_params = optimize(trials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\n",
    "             'num_round': 200,\n",
    "             'colsample_bytree': 0.65,\n",
    "             'eta': 0.145,\n",
    "             'gamma': 0.55,\n",
    "             'max_depth': 10,\n",
    "             'min_child_weight': 4.0,\n",
    "             'subsample': 1.0,\n",
    "             'eval_metric': 'mae',\n",
    "             'objective': 'reg:linear',\n",
    "             'nthread' : 4,\n",
    "             'silent' : 1,\n",
    "             'seed': 4444}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtrain = xgb.DMatrix(X_train_transformed, label=y_train_log)\n",
    "dvalid = xgb.DMatrix(X_test_transformed, label=y_test_log)\n",
    "boosting_model = xgb.train(params, dtrain, params['num_round'])\n",
    "predictions = boosting_model.predict(dvalid).reshape((X_test_transformed.shape[0], 1))\n",
    "score = mean_absolute_error(y_test, np.expm1(predictions))\n",
    "print('Boosting score: %f' % score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 8. Построение кривых валидации и обучения "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В большом количестве в предыдущем пункте"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 9. Прогноз для тестовой или отложенной выборки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В большом количестве в части 7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 10. Оценка модели с описанием выбранной метрики"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Приведем результаты различных моделей на тестовой выборке.\n",
    "Как уже оговаривалось ранее, в качестве метрики используем MAE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pure_ridge = Ridge(random_state=4444, alpha=0.00001) # гребневая регрессия на исходных данных\n",
    "pure_ridge.fit(X_train_transformed, y_train_log)\n",
    "pure_ridge_score = mean_absolute_error(y_test, np.expm1(pure_ridge.predict(X_test_transformed)))\n",
    "print('Ridge regression score: %f' % pure_ridge_score)\n",
    "poly_ridge = Ridge(random_state=4444, alpha=52.140083) # гребневая регрессия с полиномиальными признаками\n",
    "poly_ridge.fit(X_train_transformed_add, y_train_log)\n",
    "poly_ridge_score = mean_absolute_error(y_test, np.expm1(poly_ridge.predict(X_test_transformed_add)))\n",
    "print('Ridge regression score with poly features: %f' % poly_ridge_score)\n",
    "forest_score = mean_absolute_error(y_test, np.expm1(forest_model.predict(X_test_transformed)))\n",
    "print('Random forest score: %f' % forest_score)\n",
    "boosting_score = mean_absolute_error(y_test, np.expm1(boosting_model.predict(dvalid)))\n",
    "print('XGBoost score: %f' % boosting_score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Результаты близки к тем, что получались на кросс-валидации. Так что всё хорошо =)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 11. Выводы "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "В данном проекте рассматривались достаточно \"простые\" данные, поэтому основной упор был сделан на применение различных моделей для их анализа. С одной стороны, случайный лес без какой-либо настройки гиперпараметров показал лучший результат. С другой стороны, если потратить больше времени на оптимизацию градиентного бустинга, возможно, он сможет показать результат лучше, чем у случайного леса. Стоит, также, отметить линейную модель: после добавления полиномиальных признаков она показала очень неплохой результат (если сравнивать с моделью без дополнительных признаков =)). Зато сложность гораздо меньше. Если вдруг кому-то по жизни придется оценивать бриллианты, можете смело использовать предложенную модель случайного леса. В среднем будете терять по 275 $ с одного камушка :p\n",
    "\n",
    "Спасибо за внимание!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
